# Check the conda env running, it's important that you're running in the correct conda env
# Don't change anything here
Sys.getenv("CONDA_DEFAULT_ENV")

# Change the path to your own directory
rm(list=ls())
setwd('/resnick/groups/NolanLab/cchen/mannitol/cnmf/20250722_root_atlas_cnmf')

# Load the libraries
# Don't change anything here
suppressMessages(library(GeneOverlap))
suppressMessages(library(gprofiler2))
suppressMessages(library(tidyverse))
suppressMessages(library(ComplexHeatmap))
suppressMessages(library(circlize))
suppressMessages(library(Seurat))

# Check the loaded packages and their versions. 
# # Be sure to execute this! It's important for reproducibility and debugging
# Don't change anything here
sessionInfo()

R version 4.3.3 (2024-02-29)
Platform: x86_64-conda-linux-gnu (64-bit)
Running under: Red Hat Enterprise Linux 9.5 (Plow)

Matrix products: default
BLAS/LAPACK: /resnick/groups/NolanLab/cchen/miniconda3/envs/muscat/lib/libopenblasp-r0.3.27.so;  LAPACK version 3.12.0

locale:
 [1] LC_CTYPE=C.UTF-8           LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

time zone: America/Los_Angeles
tzcode source: system (glibc)

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] Seurat_5.2.1          SeuratObject_5.0.2    sp_2.2-0             
 [4] circlize_0.4.16       ComplexHeatmap_2.18.0 lubridate_1.9.4      
 [7] forcats_1.0.0         stringr_1.5.1         dplyr_1.1.4          
[10] purrr_1.0.2           readr_2.1.5           tidyr_1.3.1          
[13] tibble_3.2.1          ggplot2_3.5.1         tidyverse_2.0.0      
[16] gprofiler2_0.2.3      GeneOverlap_1.38.0   

loaded via a namespace (and not attached):
  [1] RColorBrewer_1.1-3     jsonlite_1.9.1         shape_1.4.6.1         
  [4] magrittr_2.0.3         spatstat.utils_3.1-2   farver_2.1.2          
  [7] GlobalOptions_0.1.2    vctrs_0.6.5            ROCR_1.0-11           
 [10] spatstat.explore_3.3-4 base64enc_0.1-3        htmltools_0.5.8.1     
 [13] sctransform_0.4.1      parallelly_1.38.0      KernSmooth_2.23-24    
 [16] htmlwidgets_1.6.4      ica_1.0-3              plyr_1.8.9            
 [19] plotly_4.10.4          zoo_1.8-13             uuid_1.2-1            
 [22] igraph_2.1.4           mime_0.12              lifecycle_1.0.4       
 [25] iterators_1.0.14       pkgconfig_2.0.3        Matrix_1.6-5          
 [28] R6_2.5.1               fastmap_1.2.0          fitdistrplus_1.2-2    
 [31] future_1.34.0          shiny_1.10.0           clue_0.3-65           
 [34] digest_0.6.37          colorspace_2.1-1       patchwork_1.3.0       
 [37] S4Vectors_0.40.2       tensor_1.5             RSpectra_0.16-2       
 [40] irlba_2.3.5.1          progressr_0.15.1       spatstat.sparse_3.1-0 
 [43] fansi_1.0.6            timechange_0.3.0       polyclip_1.10-7       
 [46] abind_1.4-5            httr_1.4.7             compiler_4.3.3        
 [49] withr_3.0.2            doParallel_1.0.17      fastDummies_1.7.5     
 [52] gplots_3.1.3.1         MASS_7.3-60.0.1        rjson_0.2.21          
 [55] gtools_3.9.5           caTools_1.18.2         tools_4.3.3           
 [58] lmtest_0.9-40          httpuv_1.6.15          future.apply_1.11.2   
 [61] goftest_1.2-3          glue_1.7.0             nlme_3.1-165          
 [64] promises_1.3.2         pbdZMQ_0.3-13          Rtsne_0.17            
 [67] cluster_2.1.6          reshape2_1.4.4         generics_0.1.3        
 [70] spatstat.data_3.1-4    gtable_0.3.5           tzdb_0.4.0            
 [73] data.table_1.15.4      hms_1.1.3              utf8_1.2.4            
 [76] spatstat.geom_3.3-5    BiocGenerics_0.48.1    RcppAnnoy_0.0.22      
 [79] ggrepel_0.9.5          RANN_2.6.2             foreach_1.5.2         
 [82] pillar_1.9.0           spam_2.11-1            IRdisplay_1.1         
 [85] RcppHNSW_0.6.0         later_1.4.1            splines_4.3.3         
 [88] lattice_0.22-6         deldir_2.0-4           survival_3.8-3        
 [91] tidyselect_1.2.1       miniUI_0.1.1.1         pbapply_1.7-2         
 [94] gridExtra_2.3          IRanges_2.36.0         scattermore_1.2       
 [97] stats4_4.3.3           matrixStats_1.3.0      stringi_1.8.4         
[100] lazyeval_0.2.2         evaluate_1.0.3         codetools_0.2-20      
[103] cli_3.6.3              uwot_0.1.16            IRkernel_1.3.2        
[106] xtable_1.8-4           reticulate_1.41.0.1    repr_1.1.7            
[109] munsell_0.5.1          Rcpp_1.0.13            spatstat.random_3.3-2 
[112] globals_0.16.3         png_0.1-8              spatstat.univar_3.1-2 
[115] parallel_4.3.3         dotCall64_1.2          bitops_1.0-8          
[118] listenv_0.9.1          viridisLite_0.4.2      scales_1.3.0          
[121] ggridges_0.5.6         crayon_1.5.3           GetoptLong_1.0.5      
[124] rlang_1.1.4            cowplot_1.1.3

# Read in GEP score
GEP_score_raw <- read.table("20250722_root_atlas_cnmf_parallel/20250722_root_atlas_cnmf_parallel.gene_spectra_score.k_30.dt_0_1.txt", 
                   header = T)

head(GEP_score_raw, 3)

GEP_score <- as.data.frame(t(GEP_score_raw))

head(GEP_score, 3)

colnames(GEP_score) <- paste0("GEP_", colnames(GEP_score))

head(GEP_score, 3)

GEP_score <- GEP_score %>%
    rownames_to_column("gene_ID")

head(GEP_score, 3)

sum(is.na(GEP_score))

GEP_list <- list()
for (i in seq(1,30)){
    ls <- GEP_score %>% arrange(desc(GEP_score[[i+1]])) %>%  slice_head(n = 9) %>% select(gene_ID) 
    GEP_list <- append(GEP_list, list(ls$gene_ID))
}

names(GEP_list) <- paste0("GEP_",seq(1,30))

GEP_list

## Top 9 genes for each GEP will be unique
sum(sort(table(as.character(unlist(GEP_list))))>1)

length(as.character(unlist(GEP_list)))

## For GO analysis, choose top 200 genes for each GEP
GEP_list <- list()
for (i in seq(1,30)){
    ls <- GEP_score %>% arrange(desc(GEP_score[[i+1]])) %>%  slice_head(n = 200) %>% select(gene_ID) 
    GEP_list <- append(GEP_list, list(ls$gene_ID))
}

names(GEP_list) <- paste0("GEP_",seq(1,30))

GEP_list

# Perform GO analysis
cluster_GO <- gost(GEP_list, organism = "athaliana", correction_method = "fdr", significant = F, multi_query = F)

cluster_GO_df <- cluster_GO[[1]] # the result data frame

cluster_GO_sig <- filter(cluster_GO_df, p_value<=0.01)

head(cluster_GO_sig, 3)

sort(unique(cluster_GO_sig$query))

# top  terms for each cluster

cluster_GO_sig %>%
  filter(source=="GO:BP", intersection_size>=4) %>%
  group_by(query) %>%
  top_n(5, wt = -p_value) %>%
  arrange(desc(p_value)) -> top_GO

GO_n <- cluster_GO_sig %>%
  filter(source=="GO:BP", intersection_size>=4) %>%
  group_by(term_id) %>%
  tally() %>%
  arrange(desc(n))


GO_n <- dplyr::rename(GO_n, "n_clusters"=n)

cluster_GO_sig_n <- left_join(cluster_GO_sig, GO_n)

# get all terms for the top ones so that all clusters have values

top_GO_all <- filter(cluster_GO_df, term_id %in% top_GO$term_id)


#spread and plot


top_GO_sel <- select(top_GO_all, query, p_value, term_id, term_name)

spread_GO <- spread(top_GO_sel, key = query, p_value)

spread_GO[is.na(spread_GO)] <- 1

spread_GO_m <- as.matrix(-log10(spread_GO[3:ncol(spread_GO)]))
rownames(spread_GO_m) <- spread_GO$term_name

Joining with `by = join_by(term_id)`

options(repr.plot.width = 20, repr.plot.height = 20)

GO_hm <- Heatmap(spread_GO_m, 
                 name = "-log10_pval", 
                 heatmap_legend_param = list(title_position="topcenter", color_bar = "continuous"), 
                 col = colorRamp2(c(0, 10), 
                                  c("beige", "#e31a1c")), 
                 cluster_rows = T,
                 cluster_columns = T, 
                 use_raster= FALSE, 
                 show_column_names = TRUE, 
                 show_row_names = TRUE, 
                 show_row_dend = TRUE, 
                 show_column_dend = TRUE, 
                 clustering_distance_rows = "pearson",
                 clustering_distance_columns = "pearson", 
                 row_names_gp = gpar(fontsize = 12)) 


# padding - bottom, left, top, right
draw(GO_hm, padding = unit(c(15, 15, 5, 80), "mm"), heatmap_legend_side = "left")

getwd()

head(spread_GO)

# Save the result
#write.table(spread_GO,
#            "20250722_root_atlas_cnmf_parallel.gene_spectra_score.k_30.dt_0_1_GO.txt", sep = "\t",
#            row.names = F, col.names = T, quote = F)

# The function to read Seurat object
# This is specifically for reading Seurat 3/4 objects in a Seurat 5 env
# Don't change anything here
read_seu <- function(dir) {
  seu <- readRDS(dir)
  seu@assays$spliced_RNA <- NULL
  seu@assays$spliced_SCT <- NULL
  seu@assays$unspliced_RNA <- NULL
  seu@assays$unspliced_SCT <- NULL
  seu <- suppressWarnings(suppressMessages(UpdateSlots(seu)))
  seu <- suppressWarnings(suppressMessages(UpdateSeuratObject(seu)))
  return(seu)
}

rc.integrated <- read_seu("/resnick/groups/NolanLab/shared_resources/annotated_objects/rds/integrated/Root_Atlas_v2_seu3_20220106.rds")

rc.integrated

An object of class Seurat 
71468 features across 110427 samples within 3 assays 
Active assay: integrated (17513 features, 17513 variable features)
 2 layers present: data, scale.data
 2 other assays present: RNA, SCT
 5 dimensional reductions calculated: pca, umap, umap_3D, umap_2D, umap_50

GEP_usage <- read.csv("/resnick/groups/NolanLab/cchen/mannitol/cnmf/20250722_root_atlas_cnmf/20250722_root_atlas_cnmf_parallel/20250722_root_atlas_cnmf_parallel.usages.k_30.dt_0_1.consensus.txt",
                      sep="\t")

head(GEP_usage)

# Check if the order of cell barcodes are the same
GEP_usage$X[1:5]
colnames(rc.integrated)[1:5]

# Check if the order of cell barcodes are the same
identical(GEP_usage$X, colnames(rc.integrated))

# Update metadata
rc.integrated$GEP_1 <- GEP_usage$X1
rc.integrated$GEP_2 <- GEP_usage$X2
rc.integrated$GEP_3 <- GEP_usage$X3
rc.integrated$GEP_4 <- GEP_usage$X4
rc.integrated$GEP_5 <- GEP_usage$X5
rc.integrated$GEP_6 <- GEP_usage$X6
rc.integrated$GEP_7 <- GEP_usage$X7
rc.integrated$GEP_8 <- GEP_usage$X8
rc.integrated$GEP_9 <- GEP_usage$X9
rc.integrated$GEP_10 <- GEP_usage$X10
rc.integrated$GEP_11 <- GEP_usage$X11
rc.integrated$GEP_12 <- GEP_usage$X12
rc.integrated$GEP_13 <- GEP_usage$X13
rc.integrated$GEP_14 <- GEP_usage$X14
rc.integrated$GEP_15 <- GEP_usage$X15
rc.integrated$GEP_16 <- GEP_usage$X16
rc.integrated$GEP_17 <- GEP_usage$X17
rc.integrated$GEP_18 <- GEP_usage$X18
rc.integrated$GEP_19 <- GEP_usage$X19
rc.integrated$GEP_20 <- GEP_usage$X20
rc.integrated$GEP_21 <- GEP_usage$X21
rc.integrated$GEP_22 <- GEP_usage$X22
rc.integrated$GEP_23 <- GEP_usage$X23
rc.integrated$GEP_24 <- GEP_usage$X24
rc.integrated$GEP_25 <- GEP_usage$X25
rc.integrated$GEP_26 <- GEP_usage$X26
rc.integrated$GEP_27 <- GEP_usage$X27
rc.integrated$GEP_28 <- GEP_usage$X28
rc.integrated$GEP_29 <- GEP_usage$X29
rc.integrated$GEP_30 <- GEP_usage$X30

# Chage the col name of GEP_usage
colnames(GEP_usage) <- c("cell_barcode", paste0("GEP_",seq(1,30)))

head(GEP_usage)

rc.integrated$time.anno.Li.crude <- gsub("_.*$","",rc.integrated$time.celltype.anno.Li.crude)
rc.integrated$celltype.anno.Li.crude <- gsub("Distal Columella","Distal Columella_Columella",rc.integrated$time.celltype.anno.Li.crude)
rc.integrated$celltype.anno.Li.crude <- gsub("Distal Lateral Root Cap","Distal Lateral Root Cap_Lateral Root Cap",rc.integrated$celltype.anno.Li.crude)
rc.integrated$celltype.anno.Li.crude <- gsub("Proximal Columella","Proximal Columella_Columella",rc.integrated$celltype.anno.Li.crude)
rc.integrated$celltype.anno.Li.crude <- gsub("Proximal Lateral Root Cap","Proximal Lateral Root Cap_Lateral Root Cap",rc.integrated$celltype.anno.Li.crude)
rc.integrated$celltype.anno.Li.crude <- gsub("^.*_","",rc.integrated$celltype.anno.Li.crude)
rc.integrated$time.anno.Li <- gsub("_.*$","",rc.integrated$time.celltype.anno.Li)
rc.integrated$celltype.anno.Li <- gsub("Distal Columella","Distal Columella_Columella",rc.integrated$time.celltype.anno.Li)
rc.integrated$celltype.anno.Li <- gsub("Distal Lateral Root Cap","Distal Lateral Root Cap_Lateral Root Cap",rc.integrated$celltype.anno.Li)
rc.integrated$celltype.anno.Li <- gsub("Proximal Columella","Proximal Columella_Columella",rc.integrated$celltype.anno.Li)
rc.integrated$celltype.anno.Li <- gsub("Proximal Lateral Root Cap","Proximal Lateral Root Cap_Lateral Root Cap",rc.integrated$celltype.anno.Li)
rc.integrated$celltype.anno.Li <- gsub("^.*_","",rc.integrated$celltype.anno.Li)

plot_anno <- function(rc.integrated){
order <- c("Quiescent Center", "Ground Tissue","Columella", "Lateral Root Cap", "Atrichoblast", "Trichoblast", "Cortex", "Endodermis", "Metaphloem & Companion Cell","Protophloem", "Xylem", "Procambium","Xylem Pole Pericycle","Phloem Pole Pericycle", "Protoxylem", "Metaxylem", "Unknown")
palette <- c("#9400d3", "#DCD0FF","#5ab953", "#bfef45", "#008080", "#21B6A8", "#82b6ff", "#0000FF","#e6194b", "#dd77ec", "#9a6324", "#ffe119", "#ff9900", "#ffd4e3", "#9a6324", "#ddaa6f", "#EEEEEE")
rc.integrated$celltype.anno.Li <- factor(rc.integrated$celltype.anno.Li, levels = order[sort(match(unique(rc.integrated$celltype.anno.Li),order))])
color <- palette[sort(match(unique(rc.integrated$celltype.anno.Li),order))]
DimPlot(rc.integrated, reduction = "umap", group.by = "celltype.anno.Li", cols=color, pt.size=0.2,raster=FALSE)}

options(repr.plot.width=8, repr.plot.height=6)
plot_anno(rc.integrated)

options(repr.plot.width=8, repr.plot.height=6)
DimPlot(rc.integrated, reduction = "umap", group.by = "time.anno.Li",order = c("Distal Columella","Proximal Columella","Distal Lateral Root Cap","Proximal Lateral Root Cap","Maturation","Elongation", "Transition Domain", "Proliferation Domain"),
        cols = c('#F7E7B0','#FFC400','#2B871F','#005E3B', "#deebf7", "#3182bd", '#fee0d2','#de2d26'),raster=FALSE)

Idents(rc.integrated) <- rc.integrated$time.celltype.anno.Li

plot_box <- function(annotation, value, title){
    # Prepare data
data <- data.frame(
  annotation = annotation,
  value = value
)

#Prepare mean values for each annotation and reorder
data <- data %>%
  group_by(annotation) %>%
  mutate(mean_value = mean(value)) %>%
  ungroup() %>%
  mutate(annotation = reorder(annotation, desc(mean_value)))
    
  options(repr.plot.width=24, repr.plot.height=12)
  # Create the boxplot
  print(ggplot(data, aes(x = annotation, y = value)) +
  geom_boxplot() +
  labs(x = "", y = "Usage", title = title) +
  theme_minimal(base_size=18) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) )
    
}

for (i in seq(1,30)){
    plot_box(rc.integrated$time.celltype.anno.Li, GEP_usage[,which(colnames(GEP_usage)==paste0("GEP_",i))], paste0("GEP_",i))
}

plot_feature <- function(feature){
    options(repr.plot.width=6, repr.plot.height=6)
    suppressMessages(print(FeaturePlot(rc.integrated, features = paste0("GEP_",i),pt.size=1.2,raster=FALSE, order = T)+
      scale_colour_gradient(low = "#eeeeee", high = "#ff0000")))
}

for (i in seq(1,30)){
    plot_feature(i)
}

meta_df <- rc.integrated@meta.data

head(meta_df)

sel_cols <- colnames(meta_df)[str_detect(colnames(meta_df), pattern = "GEP")]

plot_df <- meta_df %>%
  select(all_of(c("time.celltype.anno.Li", sel_cols))) %>%
  pivot_longer(
    cols = all_of(sel_cols),
    names_to = "GEP",
    values_to = "score"
  )

head(plot_df)

# Average score per GEP × time.celltype.anno.Li
# complete: Complete a data frame with missing combinations of data
# Since some cell types are missing in some time points, fill in missing combinations with score = 0
agg_df <- plot_df %>%
  group_by(GEP, time.celltype.anno.Li) %>%
  summarize(score = mean(score, na.rm = TRUE), .groups = "drop") %>%
  tidyr::complete(GEP, time.celltype.anno.Li, fill = list(score = 0))  # Fill in missing combinations with 0

head(agg_df)

agg_df$GEP <- factor(agg_df$GEP, levels = paste0("GEP_", 1:n_distinct(agg_df$GEP)))

levels(agg_df$GEP )

options(repr.plot.width = 15, repr.plot.height = 12)

ggplot(agg_df, aes(y = time.celltype.anno.Li, x = GEP, fill = score)) +
  geom_tile(color = "grey90", size = 0.3) +
  scale_fill_gradient2(low = "white", high = "firebrick3", name = "GEP usage") +
  theme_light(base_size = 12) +
  theme(
    strip.placement = "outside",
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5, size = 14),
    axis.text.y = element_text(size = 12),
    panel.grid.major = element_line(color = "grey90", linewidth = 0.1),
    panel.grid.minor = element_blank(),
    panel.border = element_blank(),
    strip.background = element_rect(fill = "darkgrey"),
    strip.text = element_text(color = "white", face = "bold")
  ) +
  labs(
    x = "",
    y = ""
  )

rc.integrated_tricycle <- readRDS("/resnick/groups/NolanLab/shared_resources/annotated_objects/rds/integrated/Root_Atlas_v2_seu4_20230907_tricycle.rds")

plot_box <- function(annotation, value, title){
    # Prepare data
data <- data.frame(
  annotation = annotation,
  value = value
)

#Prepare mean values for each annotation and reorder
data <- data %>%
  group_by(annotation) %>%
  mutate(mean_value = mean(value)) %>%
  ungroup() %>%
  mutate(annotation = reorder(annotation, desc(mean_value)))
    
  options(repr.plot.width=8, repr.plot.height=6)
  # Create the boxplot
  print(ggplot(data, aes(x = annotation, y = value)) +
  geom_boxplot() +
  labs(x = "", y = "Usage", title = title) +
  theme_minimal(base_size=18) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) )
    
}

for (i in seq(1,30)){
    plot_box(rc.integrated_tricycle$tricycleCCStage[which(rc.integrated_tricycle$time.anno.Li=="Proliferation Domain" |
                                                          rc.integrated_tricycle$time.anno.Li=="Proximal Lateral Root Cap" | 
                                                          rc.integrated_tricycle$time.anno.Li=="Proximal Columella")], 
             GEP_usage[,which(colnames(GEP_usage)==paste0("GEP_",i))][which(rc.integrated_tricycle$time.anno.Li=="Proliferation Domain" | 
                                                                            rc.integrated_tricycle$time.anno.Li=="Proximal Lateral Root Cap" | 
                                                                            rc.integrated_tricycle$time.anno.Li=="Proximal Columella")],
             paste0("GEP_",i))
}

# Check if the order of cell barcodes are the same
identical(GEP_usage$cell_barcode, colnames(rc.integrated_tricycle))

GEP_usage$cell_barcode[1]
colnames(rc.integrated_tricycle)[1]

# Update metadata
rc.integrated_tricycle$GEP_1 <- GEP_usage$GEP_1
rc.integrated_tricycle$GEP_2 <- GEP_usage$GEP_2
rc.integrated_tricycle$GEP_3 <- GEP_usage$GEP_3
rc.integrated_tricycle$GEP_4 <- GEP_usage$GEP_4
rc.integrated_tricycle$GEP_5 <- GEP_usage$GEP_5
rc.integrated_tricycle$GEP_6 <- GEP_usage$GEP_6
rc.integrated_tricycle$GEP_7 <- GEP_usage$GEP_7
rc.integrated_tricycle$GEP_8 <- GEP_usage$GEP_8
rc.integrated_tricycle$GEP_9 <- GEP_usage$GEP_9
rc.integrated_tricycle$GEP_10 <- GEP_usage$GEP_10
rc.integrated_tricycle$GEP_11 <- GEP_usage$GEP_11
rc.integrated_tricycle$GEP_12 <- GEP_usage$GEP_12
rc.integrated_tricycle$GEP_13 <- GEP_usage$GEP_13
rc.integrated_tricycle$GEP_14 <- GEP_usage$GEP_14
rc.integrated_tricycle$GEP_15 <- GEP_usage$GEP_15
rc.integrated_tricycle$GEP_16 <- GEP_usage$GEP_16
rc.integrated_tricycle$GEP_17 <- GEP_usage$GEP_17
rc.integrated_tricycle$GEP_18 <- GEP_usage$GEP_18
rc.integrated_tricycle$GEP_19 <- GEP_usage$GEP_19
rc.integrated_tricycle$GEP_20 <- GEP_usage$GEP_20
rc.integrated_tricycle$GEP_21 <- GEP_usage$GEP_21
rc.integrated_tricycle$GEP_22 <- GEP_usage$GEP_22
rc.integrated_tricycle$GEP_23 <- GEP_usage$GEP_23
rc.integrated_tricycle$GEP_24 <- GEP_usage$GEP_24
rc.integrated_tricycle$GEP_25 <- GEP_usage$GEP_25
rc.integrated_tricycle$GEP_26 <- GEP_usage$GEP_26
rc.integrated_tricycle$GEP_27 <- GEP_usage$GEP_27
rc.integrated_tricycle$GEP_28 <- GEP_usage$GEP_28
rc.integrated_tricycle$GEP_29 <- GEP_usage$GEP_29
rc.integrated_tricycle$GEP_30 <- GEP_usage$GEP_30

meta_df_tricycle <- rc.integrated_tricycle@meta.data

head(meta_df_tricycle)

meta_df_tricycle_fil <- meta_df_tricycle %>%
   filter(time.anno.Li %in% c("Proliferation Domain", "Proximal Lateral Root Cap", "Proximal Columella"))

head(meta_df_tricycle_fil, 3)

unique(meta_df_tricycle_fil$time.anno.Li)

sel_cols <- colnames(meta_df_tricycle_fil)[str_detect(colnames(meta_df_tricycle_fil), pattern = "GEP")]
sel_cols

plot_df_tri_fil <- meta_df_tricycle_fil %>%
  select(all_of(c("tricycleCCStage", sel_cols))) %>%
  pivot_longer(
    cols = all_of(sel_cols),
    names_to = "GEP",
    values_to = "score"
  )

head(plot_df_tri_fil)

# Average score per GEP × time.celltype.anno.Li
# complete: Complete a data frame with missing combinations of data
# Since some cell types are missing in some time points, fill in missing combinations with score = 0
agg_df_tri_fil <- plot_df_tri_fil %>%
  group_by(GEP, tricycleCCStage) %>%
  summarize(score = mean(score, na.rm = TRUE), .groups = "drop") %>%
  tidyr::complete(GEP, tricycleCCStage, fill = list(score = 0))  # Fill in missing combinations with 0

head(agg_df_tri_fil)

agg_df_tri_fil$GEP <- factor(agg_df_tri_fil$GEP, levels = paste0("GEP_", 1:n_distinct(agg_df_tri_fil$GEP)))

levels(agg_df_tri_fil$GEP )

options(repr.plot.width = 12, repr.plot.height = 4)

ggplot(agg_df_tri_fil, aes(y = tricycleCCStage, x = GEP, fill = score)) +
  geom_tile(color = "grey90", size = 0.3) +
  scale_fill_gradient2(low = "white", high = "firebrick3", name = "GEP usage") +
  theme_light(base_size = 16) +
  theme(
    strip.placement = "outside",
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
    axis.text.y = element_text(size = 12),
    panel.grid.major = element_line(color = "grey90", linewidth = 0.1),
    panel.grid.minor = element_blank(),
    panel.border = element_blank(),
    strip.background = element_rect(fill = "darkgrey"),
    strip.text = element_text(color = "white", face = "bold")
  ) +
  labs(
    x = "",
    y = ""
  )

plot_box <- function(annotation, value, title){
    # Prepare data
data <- data.frame(
  annotation = annotation,
  value = value
)

#Prepare mean values for each annotation and reorder
data <- data %>%
  group_by(annotation) %>%
  mutate(mean_value = mean(value)) %>%
  ungroup() %>%
  mutate(annotation = reorder(annotation, desc(mean_value)))
    
  options(repr.plot.width=8, repr.plot.height=6)
  # Create the boxplot
  print(ggplot(data, aes(x = annotation, y = value)) +
  geom_boxplot() +
  labs(x = "", y = "Usage", title = title) +
  theme_minimal(base_size=18) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) )
    
}

for (i in seq(1,30)){
    plot_box(rc.integrated_tricycle$branch.anno, GEP_usage[,which(colnames(GEP_usage)==paste0("GEP_",i))], paste0("GEP_",i))
}

# Read in the table sent by Che-Wei
topgene_new <- read.table("/resnick/groups/NolanLab/cchen/mannitol/cnmf/20250722_root_atlas_cnmf/20250722_root_atlas_cnmf_parallel.gene_spectra_score.k_30.dt_0_1_annotated.csv", 
                          sep = ",", header = T)

head(topgene_new)

# Select rank 1 ~ 200 genes for each GEP
topgene_new_top <- topgene_new %>%
   dplyr::select(1:4) %>%
   filter(Rank %in% c(1:200))

head(topgene_new_top)
nrow(topgene_new_top)

# Convert into a list
GEP_list <- topgene_new_top %>%
  group_by(GEP) %>%
  summarise(genes = list(Gene), .groups = "drop") %>%
  deframe() #  convert a two-column df into a named list

class(GEP_list)

GEP_list

length(GEP_list)
length(GEP_list[[1]])

# Perform GO analysis
cluster_GO <- gost(GEP_list, organism = "athaliana", correction_method = "fdr", significant = F, multi_query = F)

cluster_GO_df <- cluster_GO[[1]] # the result data frame

cluster_GO_sig <- filter(cluster_GO_df, p_value<=0.01)

head(cluster_GO_sig, 3)

sort(unique(cluster_GO_sig$query))

# top  terms for each cluster

cluster_GO_sig %>%
  filter(source=="GO:BP", intersection_size>=4) %>%
  group_by(query) %>%
  top_n(5, wt = -p_value) %>%
  arrange(desc(p_value)) -> top_GO

GO_n <- cluster_GO_sig %>%
  filter(source=="GO:BP", intersection_size>=4) %>%
  group_by(term_id) %>%
  tally() %>%
  arrange(desc(n))


GO_n <- dplyr::rename(GO_n, "n_clusters"=n)

cluster_GO_sig_n <- left_join(cluster_GO_sig, GO_n)

# get all terms for the top ones so that all clusters have values

top_GO_all <- filter(cluster_GO_df, term_id %in% top_GO$term_id)


#spread and plot


top_GO_sel <- select(top_GO_all, query, p_value, term_id, term_name)

spread_GO <- spread(top_GO_sel, key = query, p_value)

spread_GO[is.na(spread_GO)] <- 1

spread_GO_m <- as.matrix(-log10(spread_GO[3:ncol(spread_GO)]))
rownames(spread_GO_m) <- spread_GO$term_name

Joining with `by = join_by(term_id)`

options(repr.plot.width = 20, repr.plot.height = 20)

GO_hm <- Heatmap(spread_GO_m, 
                 name = "-log10_pval", 
                 heatmap_legend_param = list(title_position="topcenter", color_bar = "continuous"), 
                 col = colorRamp2(c(0, 10), 
                                  c("beige", "#e31a1c")), 
                 cluster_rows = T,
                 cluster_columns = T, 
                 use_raster= FALSE, 
                 show_column_names = TRUE, 
                 show_row_names = TRUE, 
                 show_row_dend = TRUE, 
                 show_column_dend = TRUE, 
                 clustering_distance_rows = "pearson",
                 clustering_distance_columns = "pearson", 
                 row_names_gp = gpar(fontsize = 12)) 


# padding - bottom, left, top, right
draw(GO_hm, padding = unit(c(15, 15, 5, 80), "mm"), heatmap_legend_side = "left")

getwd()

head(spread_GO)

# Save the result
#write.table(spread_GO,
#            "20250722_root_atlas_cnmf_parallel.gene_spectra_score.k_30.dt_0_1_annotated_by_CheWei_GO.txt", sep = "\t",
#            row.names = F, col.names = T, quote = F)

	AT1G01010	AT1G01020	AT1G03987	AT1G01030	AT1G01040	AT1G01050	AT1G01060	AT1G01070	AT1G04003	AT1G01080	⋯	AT2G07860	AT2G08040	AT2G39851	AT2G42835	AT4G03580	AT4G04815	AT4G08485	AT4G06380	AT4G21030	AT4G09045
	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	⋯	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
1	-0.0009181087	-3.582355e-05	-6.010676e-06	5.955127e-05	-3.277432e-04	0.0007601265	-3.389594e-04	-0.0003825735	-3.540773e-06	-0.0002075873	⋯	0	-6.079081e-05	0	0	-2.652584e-05	9.547121e-05	0	-3.885398e-05	0	9.241401e-06
2	-0.0004486164	2.682098e-04	-4.108028e-06	-1.261800e-04	-9.303128e-05	0.0012164017	5.476103e-05	-0.0001940035	-1.323077e-05	0.0001045421	⋯	0	3.285920e-04	0	0	1.495957e-05	-1.505614e-05	0	5.341547e-05	0	-3.415716e-06
3	-0.0009037018	-2.327463e-04	4.602577e-05	9.230912e-05	-1.327719e-04	-0.0007541299	-5.152976e-04	-0.0004510174	-2.586985e-05	-0.0002557728	⋯	0	3.557958e-06	0	0	7.854490e-05	-2.970420e-05	0	2.710244e-06	0	-1.305505e-05

	1	2	3	4	5	6	7	8	9	10	⋯	21	22	23	24	25	26	27	28	29	30
	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	⋯	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
AT1G01010	-9.181087e-04	-4.486164e-04	-9.037018e-04	5.322456e-05	1.878345e-05	-1.629325e-04	-5.377264e-04	-5.258016e-04	2.931742e-03	-6.148138e-05	⋯	-2.979482e-06	-6.110111e-04	-1.510016e-04	1.029368e-04	-2.267858e-04	3.507503e-04	-1.677782e-04	-1.491155e-04	-1.151191e-04	2.233542e-04
AT1G01020	-3.582355e-05	2.682098e-04	-2.327463e-04	-1.341628e-04	7.734951e-04	-1.128617e-04	1.238245e-04	-4.020132e-04	2.879971e-04	-1.836733e-04	⋯	2.454857e-04	-2.581357e-05	-1.463696e-04	1.477867e-04	1.615548e-04	3.322189e-04	-3.408717e-05	1.418923e-04	7.106826e-05	-2.479525e-05
AT1G03987	-6.010676e-06	-4.108028e-06	4.602577e-05	-1.051816e-05	-7.882518e-06	-8.413725e-06	-1.071538e-05	-1.024367e-05	-7.762036e-06	7.104151e-08	⋯	-8.789088e-06	-4.693212e-06	-4.245099e-06	-8.146777e-06	-1.398544e-06	9.708034e-07	-3.479765e-06	-6.700524e-06	-8.340389e-07	-8.712404e-07

	GEP_1	GEP_2	GEP_3	GEP_4	GEP_5	GEP_6	GEP_7	GEP_8	GEP_9	GEP_10	⋯	GEP_21	GEP_22	GEP_23	GEP_24	GEP_25	GEP_26	GEP_27	GEP_28	GEP_29	GEP_30
	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	⋯	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
AT1G01010	-9.181087e-04	-4.486164e-04	-9.037018e-04	5.322456e-05	1.878345e-05	-1.629325e-04	-5.377264e-04	-5.258016e-04	2.931742e-03	-6.148138e-05	⋯	-2.979482e-06	-6.110111e-04	-1.510016e-04	1.029368e-04	-2.267858e-04	3.507503e-04	-1.677782e-04	-1.491155e-04	-1.151191e-04	2.233542e-04
AT1G01020	-3.582355e-05	2.682098e-04	-2.327463e-04	-1.341628e-04	7.734951e-04	-1.128617e-04	1.238245e-04	-4.020132e-04	2.879971e-04	-1.836733e-04	⋯	2.454857e-04	-2.581357e-05	-1.463696e-04	1.477867e-04	1.615548e-04	3.322189e-04	-3.408717e-05	1.418923e-04	7.106826e-05	-2.479525e-05
AT1G03987	-6.010676e-06	-4.108028e-06	4.602577e-05	-1.051816e-05	-7.882518e-06	-8.413725e-06	-1.071538e-05	-1.024367e-05	-7.762036e-06	7.104151e-08	⋯	-8.789088e-06	-4.693212e-06	-4.245099e-06	-8.146777e-06	-1.398544e-06	9.708034e-07	-3.479765e-06	-6.700524e-06	-8.340389e-07	-8.712404e-07

	gene_ID	GEP_1	GEP_2	GEP_3	GEP_4	GEP_5	GEP_6	GEP_7	GEP_8	GEP_9	⋯	GEP_21	GEP_22	GEP_23	GEP_24	GEP_25	GEP_26	GEP_27	GEP_28	GEP_29	GEP_30
	<chr>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	⋯	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
1	AT1G01010	-9.181087e-04	-4.486164e-04	-9.037018e-04	5.322456e-05	1.878345e-05	-1.629325e-04	-5.377264e-04	-5.258016e-04	2.931742e-03	⋯	-2.979482e-06	-6.110111e-04	-1.510016e-04	1.029368e-04	-2.267858e-04	3.507503e-04	-1.677782e-04	-1.491155e-04	-1.151191e-04	2.233542e-04
2	AT1G01020	-3.582355e-05	2.682098e-04	-2.327463e-04	-1.341628e-04	7.734951e-04	-1.128617e-04	1.238245e-04	-4.020132e-04	2.879971e-04	⋯	2.454857e-04	-2.581357e-05	-1.463696e-04	1.477867e-04	1.615548e-04	3.322189e-04	-3.408717e-05	1.418923e-04	7.106826e-05	-2.479525e-05
3	AT1G03987	-6.010676e-06	-4.108028e-06	4.602577e-05	-1.051816e-05	-7.882518e-06	-8.413725e-06	-1.071538e-05	-1.024367e-05	-7.762036e-06	⋯	-8.789088e-06	-4.693212e-06	-4.245099e-06	-8.146777e-06	-1.398544e-06	9.708034e-07	-3.479765e-06	-6.700524e-06	-8.340389e-07	-8.712404e-07

	query	significant	p_value	term_size	query_size	intersection_size	precision	recall	term_id	source	term_name	effective_domain_size	source_order	parents
	<chr>	<lgl>	<dbl>	<int>	<int>	<int>	<dbl>	<dbl>	<chr>	<chr>	<chr>	<int>	<int>	<list>
1	GEP_1	TRUE	0.007873161	629	145	15	0.1034483	0.02384738	GO:0008610	GO:BP	lipid biosynthetic process	21925	3192	GO:00066....
2	GEP_1	TRUE	0.007873161	2940	145	38	0.2620690	0.01292517	GO:0042221	GO:BP	response to chemical	21925	9956	GO:0050896
3	GEP_1	TRUE	0.007873161	1097	145	20	0.1379310	0.01823154	GO:0006629	GO:BP	lipid metabolic process	21925	2334	GO:0044238

20250722 Root Atlas cNMF GEP GO analysis¶

GEP usage¶

Plot heatmaps¶

GEPs specificity on cell cycle phases¶

Plot heatmaps¶

GEPs specificity on branches¶

Redo GO analysis based on Che-Wei's annotated gene score¶

A data.frame: 6 × 32
	term_id	term_name	GEP_1	GEP_10	GEP_11	GEP_12	GEP_13	GEP_14	GEP_15	GEP_16	⋯	GEP_28	GEP_29	GEP_3	GEP_30	GEP_4	GEP_5	GEP_6	GEP_7	GEP_8	GEP_9
	<chr>	<chr>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	⋯	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
1	GO:0000278	mitotic cell cycle	1.000000000	1.00000000	1.00000000	1.000000e+00	1.0000000000	1.000000e+00	1.0000000	1.00000000	⋯	1.656482e-32	0.9860805	1.000000e+00	7.471481e-01	0.36567404	1.000000e+00	0.0103292	1.00000000	1.000000e+00	1.0000000000
2	GO:0001666	response to hypoxia	0.008568784	0.67627425	0.90535798	3.256667e-05	0.4932206146	9.649602e-41	0.3203484	0.89009059	⋯	1.000000e+00	1.0000000	3.728821e-07	6.970266e-01	0.59935653	1.963899e-14	0.9322904	0.45451009	1.000000e+00	0.1112759458
3	GO:0006082	organic acid metabolic process	0.316765040	0.10836682	0.07780494	2.811062e-01	0.4525963152	2.361847e-03	0.1299891	0.03005019	⋯	1.000000e+00	1.0000000	1.652159e-03	1.687448e-08	0.38597795	1.234134e-01	0.1195624	0.01329521	4.443228e-17	0.0088613142
4	GO:0006259	DNA metabolic process	1.000000000	0.99129690	1.00000000	1.000000e+00	1.0000000000	1.000000e+00	0.9969347	1.00000000	⋯	4.673976e-01	0.9593055	9.954154e-01	1.000000e+00	0.98144213	9.805975e-01	0.9729743	0.99280585	1.000000e+00	0.9844907236
5	GO:0006260	DNA replication	1.000000000	0.78169780	1.00000000	1.000000e+00	1.0000000000	1.000000e+00	1.0000000	1.00000000	⋯	5.545803e-02	0.9860805	1.000000e+00	1.000000e+00	0.71923409	1.000000e+00	0.8036567	1.00000000	1.000000e+00	1.0000000000
6	GO:0006355	regulation of DNA-templated transcription	0.491949051	0.08296448	1.00000000	1.000000e+00	0.0001186245	1.647069e-04	0.6833902	1.00000000	⋯	1.000000e+00	1.0000000	7.128563e-01	2.481224e-01	0.08138804	1.828872e-01	0.9945004	0.19307211	9.516024e-01	0.0001418159

A data.frame: 6 × 31
	X	X1	X2	X3	X4	X5	X6	X7	X8	X9	⋯	X21	X22	X23	X24	X25	X26	X27	X28	X29	X30
	<chr>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	⋯	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
1	AAACCCAAGCCATTCA_1	0.04893713	0.090610156	0.007496949	0.004143808	0.02813839	0.03238366	0.44500112	0.003822049	0.00000000	⋯	0.000000000	0.003786700	0.005426142	0.003614915	0	0.003510636	0.0002423366	0.00000000	0.000000e+00	0.0001246741
2	AAACCCAAGCTGTTCA_1	0.00000000	0.000000000	0.000000000	0.009822821	0.57054346	0.00000000	0.10012456	0.011877956	0.00000000	⋯	0.000000000	0.000000000	0.003494375	0.003057198	0	0.000000000	0.0042506721	0.00000000	1.512234e-05	0.0000000000
3	AAACCCAAGGCGAACT_1	0.00000000	0.000000000	0.000000000	0.000000000	0.14837877	0.54437984	0.21511212	0.000000000	0.21121791	⋯	0.000000000	0.000000000	0.000000000	0.000000000	0	0.000000000	0.0000000000	0.00000000	0.000000e+00	0.0000000000
4	AAACCCACACCAGCCA_1	0.10637848	0.010682129	0.549489400	0.000000000	0.07846621	0.01009298	0.00000000	0.002775509	0.03702748	⋯	0.008774841	0.011675368	0.000000000	0.001699086	0	0.002404297	0.0038288122	0.00000000	0.000000e+00	0.0000000000
5	AAACCCACACTGGCGT_1	0.00000000	0.008458047	0.916310036	0.000000000	0.00000000	0.00000000	0.00000000	0.000000000	0.00000000	⋯	0.000000000	0.002467955	0.009840735	0.000000000	0	0.000000000	0.0016002835	0.00154616	1.519888e-03	0.0000000000
6	AAACCCACAGACCTGC_1	0.40810441	0.148901553	0.103435620	0.000000000	0.07170737	0.00000000	0.07358418	0.000000000	0.04873583	⋯	0.000000000	0.000000000	0.002350374	0.001118176	0	0.000000000	0.0009722757	0.00000000	0.000000e+00	0.0002402171

A data.frame: 6 × 31
	cell_barcode	GEP_1	GEP_2	GEP_3	GEP_4	GEP_5	GEP_6	GEP_7	GEP_8	GEP_9	⋯	GEP_21	GEP_22	GEP_23	GEP_24	GEP_25	GEP_26	GEP_27	GEP_28	GEP_29	GEP_30
	<chr>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	⋯	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
1	AAACCCAAGCCATTCA_1	0.04893713	0.090610156	0.007496949	0.004143808	0.02813839	0.03238366	0.44500112	0.003822049	0.00000000	⋯	0.000000000	0.003786700	0.005426142	0.003614915	0	0.003510636	0.0002423366	0.00000000	0.000000e+00	0.0001246741
2	AAACCCAAGCTGTTCA_1	0.00000000	0.000000000	0.000000000	0.009822821	0.57054346	0.00000000	0.10012456	0.011877956	0.00000000	⋯	0.000000000	0.000000000	0.003494375	0.003057198	0	0.000000000	0.0042506721	0.00000000	1.512234e-05	0.0000000000
3	AAACCCAAGGCGAACT_1	0.00000000	0.000000000	0.000000000	0.000000000	0.14837877	0.54437984	0.21511212	0.000000000	0.21121791	⋯	0.000000000	0.000000000	0.000000000	0.000000000	0	0.000000000	0.0000000000	0.00000000	0.000000e+00	0.0000000000
4	AAACCCACACCAGCCA_1	0.10637848	0.010682129	0.549489400	0.000000000	0.07846621	0.01009298	0.00000000	0.002775509	0.03702748	⋯	0.008774841	0.011675368	0.000000000	0.001699086	0	0.002404297	0.0038288122	0.00000000	0.000000e+00	0.0000000000
5	AAACCCACACTGGCGT_1	0.00000000	0.008458047	0.916310036	0.000000000	0.00000000	0.00000000	0.00000000	0.000000000	0.00000000	⋯	0.000000000	0.002467955	0.009840735	0.000000000	0	0.000000000	0.0016002835	0.00154616	1.519888e-03	0.0000000000
6	AAACCCACAGACCTGC_1	0.40810441	0.148901553	0.103435620	0.000000000	0.07170737	0.00000000	0.07358418	0.000000000	0.04873583	⋯	0.000000000	0.000000000	0.002350374	0.001118176	0	0.000000000	0.0009722757	0.00000000	0.000000e+00	0.0002402171

A data.frame: 6 × 135
	orig.ident	nCount_RNA	nFeature_RNA	nCount_spliced_RNA	nFeature_spliced_RNA	nCount_unspliced_RNA	nFeature_unspliced_RNA	percent.mt	percent.cp	nCount_SCT	⋯	GEP_23	GEP_24	GEP_25	GEP_26	GEP_27	GEP_28	GEP_29	GEP_30	time.anno.Li	celltype.anno.Li
	<chr>	<dbl>	<int>	<dbl>	<int>	<dbl>	<int>	<dbl>	<dbl>	<dbl>	⋯	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<chr>	<chr>
AAACCCAAGCCATTCA_1	sc_12	11727	4244	10032	3899	1695	1027	0.2899292	0.07674597	18223	⋯	0.005426142	0.003614915	0	0.003510636	0.0002423366	0.00000000	0.000000e+00	0.0001246741	Transition Domain	Atrichoblast
AAACCCAAGCTGTTCA_1	sc_12	23585	4301	20531	4064	3054	1123	0.1823193	0.03391986	19971	⋯	0.003494375	0.003057198	0	0.000000000	0.0042506721	0.00000000	1.512234e-05	0.0000000000	Elongation	Atrichoblast
AAACCCAAGGCGAACT_1	sc_12	25828	5211	22837	4885	2991	1339	2.8147747	0.38330494	20252	⋯	0.000000000	0.000000000	0	0.000000000	0.0000000000	0.00000000	0.000000e+00	0.0000000000	Elongation	Cortex
AAACCCACACCAGCCA_1	sc_12	10172	3155	8968	2917	1204	742	1.6220999	0.06881636	18426	⋯	0.000000000	0.001699086	0	0.002404297	0.0038288122	0.00000000	0.000000e+00	0.0000000000	Distal Columella	Columella
AAACCCACACTGGCGT_1	sc_12	7883	2387	7051	2199	832	523	1.5349486	0.15222631	17877	⋯	0.009840735	0.000000000	0	0.000000000	0.0016002835	0.00154616	1.519888e-03	0.0000000000	Distal Columella	Columella
AAACCCACAGACCTGC_1	sc_12	16718	4885	14501	4566	2217	1183	1.3817442	0.11364996	18505	⋯	0.002350374	0.001118176	0	0.000000000	0.0009722757	0.00000000	0.000000e+00	0.0002402171	Proximal Lateral Root Cap	Lateral Root Cap